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Abstract 



Reciprocal space methods for solving Poisson's equation for finite charge 
distributions are investigated. Improvements to previous proposals are pre- 
O ■ sented, and their performance is compared in the context of a real-space den- 

sity functional theory code. Two basic methodologies are followed: calcula- 
tion of correction terms, and imposition of a cut-off to the Coulomb potential. 
Qh| We conclude that these methods can be safely applied to finite or aperiodic 

systems with a reasonable control of speed and accuracy. 
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I. INTRODUCTION 



Density-functional theoryEm in its time-dependentH as well as ground or time-independent 
forms has proved to be an efficient method for treating electron-electron interactions and 
has been applied successfully to finite systems such as clusters,@ to bulk systems or surfaces, 
and to aperiodic systems such as defects.! However, the high computational cost of treating 
large systems places a practical limit on the size of systems that can be studied. 

The use of pseudopotentialsi@ enhances the performance of this sort of calculation by 
avoiding an explicit treatment of the Kohn-Sham orbitals associated with the core. Further- 
more, the smoothness of the resulting valence pseudowavefunctions allows the use of a plane 
wave basis for describing them, and consequently also the electronic density. A plane wave 
basis is particularly attractive because it allows use of the fast fourier transform (FFT), for 
rapid and memory efficient transformations. 

A discrete but truncated set of plane waves based on the reciprocal lattice is one natural 
basis set for a periodic system. However, for finite systems, or more generally, for systems 
lacking periodicity such as defects in solids, the use of a discrete set of plane waves will 
generate periodic images of the finite cell to be studied. In the case of a finite system this 
leads to a problem in the calculation of the electrostatic potential due to the electrons, the 
so-called Hartree potential, due to the long range of the Coulomb interaction. Nevertheless, 
discrete plane wave basis sets are often used for finite systems because of the great efficiency 
of the FFT, and errors in the Hartree potential due to periodic images are usually ignored, 
or reduced by increasing the size of the supercell. These spurious effects might seriously 
affect the calculated equilibrium structure and dynamics of weakly bounded molecules or 
clusters, eg. water. However, several methods have been proposed recently for treating this 
problem.Biii 

Our purpose here is to compare four methods for solving Poisson's equation for finite 
systems. One of them is an iterative, real-space method based on finite differences and 
conjugate-gradients minimization, which obviously doesn't suffer from the problems related 
to periodic images. The other three use discrete plane wave basis sets and FFT's, but treat 
the cell-to-cell interactions in different ways. Two of these plane wave methods impose 
a cut-off to the Coulomb interaction in real space, and have been described elsewhere .0 
However, for one of these, which uses what we term a cubic cut-off, we have found significant 
improvements which we believe will be of interest to practitioners. The third plane wave 
method has been developed and tested by us, but we have found a close relationship between 
it and the local moment countercharge (LMCC) method proposed by Schultz,0 and also to 
the Fourier analysis with long range forces (FALR) method of Lauritsch and Reinhard.i 
However, our scheme is formulated more generally, and allows for better control of errors. 

In order to compare the performance of the different methods we have studied some 
exactly soluble model systems, and NaCl and Na^ 2 molecules which, because of their polar 
or charged nature, are difficult to treat with plane wave methods. Of particular interest 
is the way the computational time scales with system size. Although all the plane wave 
methods scale as a few times N log N, where N is the number of real-space mesh points, the 
proportionality factor varies substantially from method to method. We shall compare the 
speed and memory requirements of the methods, and how these scale with the size of the 
systems. 



2 



The plane wave schemes we discuss are intended mainly to deal with neutral or charged 
molecules or clusters in free space and could be implemented directly in existing ab initio 
plane wave or real-space codes, but they could also be used in those LCAO basis set codes 
which base the calculation of the Hartree potential on the FFT. General subroutines for 
calculating the Hartree potential using these methods are available upon request from the 
authors or can be downloaded from the web page.0 

Short theoretical descriptions of the plane wave methods are given in section II where 
we emphasize the improvements we have developed. Section III presents and compares 
the results for the Hartree potential calculated using the different methods, and concluding 
remarks are made in the final section. Atomic units are used unless otherwise stated. 



II. THEORETICAL DESCRIPTION OF PLANE-WAVE METHODS 

A. Uncorrected calculations 

The solution of Poisson's equation, V 2 V^ + Airn = 0, which goes to zero at infinity, for 
a charge density n, localized within a cell C of volume Q, is given by 

V H [n,r]= [ dv'^L. (1) 
J |r -r| 

Within the cell n and Vh may be expressed as Fourier series: n (r) = Cr x ^y Gr n(G) 

G 

where n(G) = J drn(r)e~ tGr , and similarly for Vjj, and where the G vectors are reciprocal 

vectors of the lattice formed by repeating the cell C. If the Fourier coefficients, n(G), are 
negligible for G larger than some cut-off so that the sums over G may be truncated, then the 
n(G) and the n(r) points are related through a discrete Fourier transform. This amounts 
to approximating the integral over the cell in the definition of n(G) by the trapezium rule, 
a point to which we shall return later. However, n(r) given by the Fourier sum is periodic 
so that the straightforward substitution into Eq. ([!]) gives a potential: 

n».r] = |£^ & (2) 

G^O 

which differs from Vh- But the merit of V is that it can be calculated using the the very 
efficient FFT with its NlogN scaling, and so we now modify or apply corrections to Eq. 
so that it can be used to obtain Vh- Two aspects of V given by Eq. (||) require attention. 

• The G = component in Eq. @ is arbitrarily set to zero. For a charged system, this 
corresponds physically to introducing a uniform compensating charge background, b, 
so that the system is electrically neutral. For a neutral system, this means that the 
boundary condition, V(r — > oo) = 0, is not satisfied (which, of course, also happens 
in the charged case). 

• V is the potential due to the charge distribution n in the central cell plus that due to 
the images of n — b in all other cells. 
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Since we are dealing with the electron charge distribution, there is obviously a net charge. 
The fact that the whole system (cores + electrons) may or may not be charged, is irrelevant 
for the discussion presented in this paper. However, it might be important if the calculations 
are total-energy supercell calculations. In this case, the system of ion cores is also treated 
using reciprocal space, so that a background of opposite sign has to be added. If the 
finite system is neutral, the effect of the backgrounds cancels. The spurious effect of higher 
multipoles, however, will remain. The distinction is important, though, because the uniform 
background introduces an important error in the total energy of order 0(L _1 ), L being the 
size of the cell, whereas the leading effect of the presence of the multipoles is the dipole-dipole 
term, behaving like 0(L~ 3 ). This is shown in the calculations presented in next section. 



We can start to deal with the cell-to-cell interaction, by eliminating the effects of net 
charge. This can be done by subtracting from the original charge distribution, n, an auxiliary 
charge distribution n aux , so that no net charge remains. The potential Vh then becomes: 



The term Vn[n — n aux \ can be treated using the FFT techniques, and then the correction 
Vff[n aux ], calculated explicitly in real space, added on. This method is especially convenient 
if the Fourier components of n aux can be calculated analytically, so that the expression 
becomes: 



where V[n] follows the definition given in Eq. (Q), and ip = V[n aux ] — Vu\n aU x\ is a function 
which can be calculated analytically. The sign rs denotes that the effect of higher multipoles 
is still included. The choice of n aux is arbitrary; it could be a uniform density, or a eaussian 
centered on the origin, in both cases the function ip can be calculated analytically.E3 

This procedure can be easily generalized, to account for the interaction of higher multi- 
poles. We merely need to add an auxiliary charge distribution which mimics the multipoles 
whose effect we want to subtract. This procedure is called by Shultzll local moment counter- 
charge method (LMCC). Shultz accounts for the monopole and dipole corrections through 
a superposition of localized Gaussian charge distributions constructed to have the same net 
charge and dipole moment. Higher multipoles can similarly be accounted for by the super- 
position of additional Gaussian distributions, but the procedure becomes complicated. A 
straightforward approach which is more easily generalized introduces an auxiliary charge 
distribution in the form: 



B. Multipoles-corrections method 



V H [n] = V H [n - n aux ] + V H [n aux \. 



(3) 



V H [n]*V[n]-ij), 



(4) 



oo 




(5) 



Z=0 m=-l 



where 
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and Mi m is the multipole moment of n given by 

M lm = [ drn{r)r l Y lm {v). (7) 



The width parameter a is to be chosen so that n aux is negligible at the cell boundary. If 
high order moments are required a could be taken to be /-dependent, decreasing somewhat 
with I. Note that the I = term in Eq. (|5|) corrects for the net charge as described above. 

The auxiliary density is localized within a cell and has the same multipole moments as n. 
We can now correct for the presence of the periodic images of n, and obtain for the required 
Hartree potential in the central cell: 

oo I 

V H [n, r] = V [n, r] - ^ ^ M lm ^ lm + V , (8) 

1=0 m=-l 

where Vq is a constant shift yet to be determined, and the functions ipim, which are inde- 
pendent of the charge distribution n, are given by 

= ^pFTT^ g G'-V" 3c2/4 > U G)e— _ V^_i_ /l(r/o)1Ur) . (9) 

The first term in ipi m is the periodic potential due to ni m in every cell, and the second 
term subtracts the effect of ni m in the central cell. The function l\ is: 

I x {x) = ^ dtt 2l e~ t2 . (10) 



.7 

The procedure for obtaining V# is to calculate and store the functions ipi m once and for 
all for as many of the multipoles as are needed to achieve the desired precision, then for the 
particular charge distribution V is calculated using FFTs, the Mi m are computed, and Eqs. 
(§) and (|9]) used to obtain Vh within the central cell apart from the constant shift. 

The shift Vq is chosen so that the boundary condition V^(r — > oo) = is satisfied. This 
we accomplish by computing the average value over the surface of the cell of the corrected 
but unshifted potential Vh obtained by simply putting Vq = in Eq. (H). For a cubic cell 
we have:0 

Since n is zero at the boundary we may interchange the order of integration in Eq. (|TT|) 
to give: 

Vh = J dr'n(r')u(r'), (12) 



where 



u(r') 



ds 1 (13) 



c 
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is the potential inside the cube due to a unit charge uniformly distributed over the cube 
surface. As such, by Gauss's theorem and symmetry, u is constant inside the cell and has 
the value a/Vt 1 ^ with a = 1.586718 for a cube. If n integrates to z, then we have finally: 

W=^- (14) 

and Vq should be chosen so that Vh has this value. The calculation of the correct shift to 
satisfy the boundary condition requires the computation of the surface average of V, but 
since the Kohn-Sham orbitals are unaffected by the addition of a constant to the potential, 
this only needs to be performed at the end of a self-consistent cycle.0 



C. Cut-off methods 



We now review two established methodsp" based on imposing a cut-off on the Coulomb 
interactions. They are exact, but need a bigger cell, which is a computational drawback as 
we shall see. 

Let us define a new cell D, which includes C. This new cell will define new coefficients 
Go, which are the reciprocal vectors of the lattice formed by repeating D. Retrieving Eq. 
(|), we realize that the function V[n] can also be expressed as: 

V\n,r]= [ dr'n {p \r')—^—, (15) 
J |r' — r| 

where is the function formed by the sum of n and all its periodic repetitions in the 
superlattice. 

We now introduce a truncated Coulomb potential, with the following properties: 

.. f , ,, , for r and r' both belonging to the same image of C. , 
/(r — r ) = < i rr i (16) 
^ , for r and r' belonging to different images of C. 

It is easily seen that: 

J dr'n ip) {r')f(r' - r) = V H [n, r] (17) 
for every r G C. And thus we can calculate Vh as we did for V in Eq. (^): 

4-7T 

Vn(r) = -Y,HG)UG)e lGr , (18) 

G^O 

where /(G) = / drf(r)e~ iGr . 

Two choices for D and / have been given. One of them uses a spherical shapei for the 
cut-off of the Coulombic interaction, and the other a cubic shape,0 based on the assumption 
that the original cell C is cubic itself. 
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1. Spherical cut-off method 



Let L c be the length of the side of C. We will define a larger cubic cell of side L D = 
(1 + \/3)Lc, centered on C. For this choice of D we define next the truncated Coulomb 
interaction: 



Defined in this way, / meets the required conditions expressed in Eq. fllq) , because any 
two points belonging to C are always closer than \/3Lc, and any two points belonging to 
different images of C are always farther away than y/3Lc- 

The Fourier transform of n has to be calculated numerically in the larger cell D; however 
that of / is easily obtained analytically: 

m(G )=4.I^*lM. (20 ) 



2. Cubic cut-off method 

The former proposal is exact; but a very large cell is needed, which increases the time to 
evaluate the FFTs. Reducing Lc introduces spurious interactions and thus spoils the pre- 
cision of the calculations, but if extremely precise calculation are not needed, a compromise 
could be reached. 

Our aim now is to reduce Lp but maintain an accurate evaluation of Vjj. We take the 
larger cell to have Ld = 2Lc and the cut- off Coulomb interaction to be: 

f(r - r') = ( ' r 7 r ' G ° (21) 

[ , otherwise. 

If r and r' belong to C, r — r' belongs to D. And if r and r' belong to different images 
of C, then r — r' will not belong to D. Thus again / is correctly defined. 

The Fourier transform of this function / has to be calculated numerically, and we face 
here two drawbacks: the function has a singularity at the origin, and is not analytic at the 
boundary. 

Jarvis et afl dealt with the first problem by integrating and averaging the singularity 
over a grid unit, which may not be adequate. The second problem appears to have been 
overlooked. In any event, their treatment of the Mg atom using the cubic cut-off method 
converges poorly compared with the results with the spherical cut-off, and they declare a 
preference for the spherical cut-off despite the larger cell size required. However, we shall 
show how to overcome these difficulties so that the cubic cutoff method, to be preferred 
because of its smaller cell size, can be used with great precision. 

The singularity at r = r' can be treated as follows: 

, 1 „-iGr_ f j_erf(r/a) ^ iGr , f ^ 1 - erf (r/a) ^_ iGr 



dr-e-^ r = / dr v / '- e~^ r + / dr L^ e " iUr , (22) 

d r Jd r J D r 
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where a is chosen small enough so as to make 1 — erf (r/a) negligible at the cell boundaries. 
The second term can be calculated analytically 

/ dr l - e f r/a) e-^ = J{1 - e"^} (23) 

and the numerical integration reduces to the first term, which is free of singularities. Even 
so, this term cannot be calculated by simply applying an FFT because the repeated function, 
although periodic, is not analytic at the boundary. Use of the FFT amounts to using the 
trapezium rule for the integration, which is exact for a periodic analytic function, but leads 
to substantial errors when there are discontinuous derivatives as we have in this case. We 
evaluated the integral using a second-order Filon's method,^! which proved to be effective. 
Other procedures^ (Simpson's, Romberg's...) could have been used - they are all rather slow 
if accurate results are to be obtained, but this calculation needs to be done only once for a 
cubic cell. If we denote by I[D(L)] (m, n 2 , n.3) the integral in Eq. ( |2"2"| ) for a cubic box of side L 
and frequency indices (%, n 2 , n.3), it is clear that I[D(L)](ni, n 2 , n 3 ) = L 2 J[D(l)](n 1 , n 2 , n 3 ). 

III. RESULTS 

A. Exactly soluble systems 

It is interesting to see the effect of the various multipoles that a charge distribution 
might have by using the multipoles-correction method on an exactly soluble system. We 
have studied systems consisting of superpositions of Gaussian charge distributions placed at 
various points Rj within the cubic cell of side L: 

exp(— J — P-) 

"( r ) = E* ^3/2 • ( 24 ) 

We have investigated the efficiency with which the multipole corrections remove the 
effects of the images of the charge distributions in other cells. This has been done as 
functions of L, as for a large enough cell the results for the periodic system should become 
exact, but at rates depending on the order of the multipoles. The results are shown in Fig. [I] 
for the cases in which there are (i) no corrections (by which we mean that only the constant 
to meet the proper boundary conditions is added to the raw potential obtained from the 
Fourier transform), (ii) monopole corrections, (iii) monopole + dipole corrections, and (iv) 
monopole + dipole + quadrupole corrections. The following points are noteworthy. 

• There is a serious, roughly 10%, error in the total energy when the Hartee potential is 
uncorrected. Although this is not a consideration in superlattice calculations provided 
the system is neutral, it is an important matter in real space calculations when the 
Hartree potentials due to the electrons alone is calculated in reciprocal space. 

• The time for the calculations behaves roughly as 0(L 3 logL), but with irregularities. 
The efficiency of the FFT algorithm depends on the prime factorization of the number 
of points to be transformed. The original FFT was developed for powers of two, but 
now algorithms exist with more flexibility.§§ We have used the FFTW package,Ezl 
with support for all the primes involved in our calculations. 
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• Adding the quadrupole corrections does not seem to improve the accuracy of the result, 
nor is the L-dependence improved. This is because the interaction energy between the 
dipole of the charge distribution in the central cell and octupoles in other cells, has 
the same L-dependence and order of magnitude as the quadrupole-quadrupole energy. 
Consequently, although the potential will be improved by adding to it the quadrupole 
corrections, there could be no significant improvement in the total energy if the system 
has a strong dipole. In general, it can be shown that the error in the electrostatic energy 
due to the presence of an I multipole in the charge distribution in the central cell, and 
I' multipoles in all other cells goes to zero like L~^ l+l +1 \ or in some special cases faster 
due to symmetry (for instance, if I = and I' < 3, or I and /' have different parity). 
Thus, adding octupole corrections to the potential will not change the L-dependence 
of the total energy if the system is charged because of the interaction of the monopole 
with the / = 4 multipole. Our calculations below on the Na^~ 2 cluster provide an 
interesting example of this behaviour. 



B. Real systems 

We have performed several electronic structure calculations on real systems to assess the 
performance of the methods. We have used a real-space codejEI in which a superlattice and 
plane waves are only used to accelerate the solution of Poisson's equation for the electron 
charge distribution. In this type of approach a correction for the net charge is always needed 
irregardless of whether the molecule or cluster itself is charged or neutral. Furthermore, in 
this approach the value of the multipoles will depend on the position of the molecule with 
respect to the centre of the cell. In order to minimize the multipole corrections the centre 
of charge of the system of ions should be placed at the centre of the cell. If this is not done 
in real space calculations the errors caused by cell-to-cell interactions could be magnified. 
In order to illustrate the effects we take the center of charge as the cell center for one of our 
test cases, and not for the other. As for other details of the calculations we used density- 
functional theory with the local-density approximation for exchange and correlation, and 
Troullier-Martinsc3 nonlocal, norm-conserving pseudopotentials. 

Our first choice for a realistic system was the NaCl molecule, also treated by ShultzS 
and Jarvis et a/,EJ because of its strong dipole moment (experimental value of 9.0D in the 
gas phase, as reported by Nelson et a/.E3) In this case, the center of charge of the system of 
ions is placed at the center of the cell. The equilibrium bond-length was calculated: (i) using 
the spherical cut-off method which is exact with a large enough cutoff, and (ii) using the 
multipoles correction and correcting only for the monopole term so as to show the influence 
of the dipole-dipole interactions which are ignored. Our calculated "exact" value is 2.413A, 
whereas the result ignoring the dipole-dipole interactions is 2.448A. 

Next, we investigated the performance and accuracy of the methods by determining 
the errors in the total energy and electric dipole moment, while monitoring the calculation 
times. We compared results for the energy and dipole moment against those obtained with 
the spherical cutoff method with a cut-off radius of y/3L c , grid parameter (0.2A) and cell 
size (L=10A). In this way an electric dipole moment of 8. 455 ID was obtained. Each of 
the four methods was then used to converge the electronic ground state of the molecule for 
successive values of a "control parameter" for speed and accuracy: 
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• For the real-space, conjugate-gradients method this parameter was the order of the 
difference formula used to evaluate derivatives. 

• For the spherical cut-off method, we note that, if the electron density is well localized 
within the C cell, the need for the full cut-off radius, \/3Lc, may be relaxed and 
a correspondingly smaller D cell used, introducing some error but accrueing time 
savings. We have investigated the effect of using a reduced cut-off radius, r cut _ fr, 
through a control parameter, a, which is the ratio of the D and C cubic cell edges: 

Ld -. . ^cut-off /or\ 

a = — = l + — (25) 

L/C L,c 



so that a = 1 + \/3 is the minimum value for which exact results are guaranteed. 

• The D cell size can also be reduced in the case of the cubic cut-off method, and the 
control parameter is again a = where a > 2 guarantees exact results. 

• For the multipoles correction method, the order of the multipoles corrected for is the 
control parameter. 

In Fig. ^ we illustrate the results obtained for each of the methods. Both cut-off methods 
are presented in the same column as they use the same control parameter, although the 
ranges of values are different. 

1. The real-space method is significantly slower than the other methods for the same 
accuracy, and a case can be made for using reciprocal space methods for calculating 
the Hartree potential in what are otherwise real-space codes. However, enhancements 
of the conjugate-gradients method are possible through preconditioning and multigrid 
techniques.@0 

2. The cut-off methods reach acceptable accuracy much below the values for a which 
guarantee exact results: 1 + \/3 and 2 respectively, for the spherical and cubic cut- 
off. This to be expected when the charge distribution is well localized within the cell. 
However, it is clearly demonstrated that, for a given accuracy, the size of the auxiliary 
cell is smaller for the cubic cut-off method, and as a result, the calculation time is also 
shorter. 

3. The multipoles-correction method already gives good accuracy if the dipole interac- 
tions are corrected for (5 x 10~ 5 eV error in the energy, and 10 _4 D error in the electric 
dipole). Without the dipole correction, the error in the energy is 0.085eV, and in 
the dipole is 0.17D, which give an indication of the size of errors to be expected when 
supercell calculations are performed for neutral molecules and no corrections are made. 

We have also performed calculations on the Na+ 2 cluster containing the same number 
of valence electrons as the NaCl molecule. Results are similar to those presented for NaCl, 
but some differences should be reported. In this case the center of charge was not placed 
at the center of the cell, consequently, although the cluster has a calculated net dipole of 
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4.5D, the electronic dipole responsible for the errors is a much larger 10. 2D. The cluster 
was positioned in the cell so that the charge density occupied most of the cell, allowing an 
optimally small cell. As a result, to achieve similar accuracy as for the NaCl molecule, we 
should expect the need for (i) larger cut-off lengths for the cut-off methods, and (ii) higher 
multipole corrections. 

In Fig. § we show the error in the energy obtained by using the multipoles correction 
method. It can be seen how the inclusion of the dipole correction yields a much less satis- 
factory error in the energy than for NaCl. Furthermore, for the reasons given earlier, the 
removal of the quadrupole-quadrupole, dipole-octupole, and octupole-octupole terms does 
not significantly improve the accuracy. Only by including all corrections to the potential up 
to fourth order multipole moments do we obtain a comparable result for the energy. The 
calculation time, which is also shown in the figure, is beginning to increase sharply by fourth 
order as further corrections are added. 

In Fig. |] we present, as well, the results for the error in the total energy and the calcu- 
lation time for the cubic-cutoff method. Comparison with Fig. [| confirms that the energy 
converges much less rapidly as a function of a than for the NaCl molecule. 



IV. CONCLUSIONS. 

We have studied some of the methods which have been proposed recently for solving 
Poisson's equation in reciprocal space for electronic structure calculations on finite systems. 
We also propose a method based on multipole corrections. Test calculations have been 
performed to assess the performance of the methods. We conclude that reciprocal-space 
methods can be accurate enough for finite or aperiodic systems, and their efficiency is a 
significant improvement over that of real-space methods. Two basic reciprocal-space meth- 
ods have been investigated: one which imposes a cut-off on the Coulomb potential, and 
one based on the removal of the spurious effects through a multipole expansion. Both yield 
satisfactory results, and comparable efficiency. 

The former approach has been already been surveyed by Jarvis et al .0 There are two 
possibilities for the cut-off function, one, the spherical cut-off, was highlighted for having 
superior convergence with the plane wave cut-off of the reciprocal lattice. However, we 
have identified and corrected problems with the other possibility, the cubic cut-off, which 
eliminates the poor behaviour, and makes this cut-off method the better of the two because 
smaller FFT's are allowed. 

The method based on multipoles corrections was initiated by Shultzjll but we have 
developed a scheme which we think is more general. Its performance is more predictable 
than that of the cut-off methods, which are sensitive to the choice of the cut-off length. The 
reason for the sensitivity is that the cut-off length determines the size of a larger auxiliary 
cell and the number of grid points over which FFT calculations are performed, and the 
FFT is sensitive to the prime number decomposition of the number of points. On the other 
hand, the speed and accuracy of our "multipoles correction" method, are adequate for most 
applications, and can be easily controlled by choosing the order of corrections applied. 

All the methods have been presented assuming a cubic cell. However, generalizations 
to other cell shapes are possible if the geometry of the system requires it. The multipoles 
correction method is immediately applicable to any cell. Clearly the spherical cut-off method 
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would be inefficient for elongated cells because the radius of the cut-off sphere is determined 
by the longest dimension of the cell. But, the cubic cut-off method can easily be generalized 
to other cell shapes, at the cost of more, and more lengthy calculations of the Fourier 
transforms of the truncated Coulomb interaction. 

We have made a simple implementation of the solvers within the self-consistent frame- 
work, but smarter algorithms can be developed, since not all the iterations of a self-consistent 
calculation need be done with the same accuracy. For example, significant improvements 
in efficiency can be gained if, for a given method, the iterations are started with a fast but 
inexact solver through appropriate choice of the control parameter, but improved as self- 
consistency is approached. Moreover, methods could be combined using, for instance, the 
real-space method for the last few iterations because of its efficiency when a good starting 
point is known. 
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FIG. 1. Error in the electrostatic energy for the system of Gaussian charges, Eq. (24), 
(continuous line) and total time of calculations (dashed line), for the indicated order used of the 
multipoles correction. 
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FIG. 2. Error in the electrostatic (first row), electric dipole (second row) and time of cal- 
culations (third row) for the NaCl molecule, using the methods indicated, as a function of the 
respective "control parameter" (see text). For the cut-off methods, crosses refer to the spherical 
cut-off method, and circles to the cubic cut-off method. All scales are logarithmic. 
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FIG. 3. Error in the electrostatic energy (continuous line) for the Na+ 2 cluster, using the 
multipoles correction method, as a function of the order of corrections included in the calculations. 
Also shown is the time of calculation for each case (dashed line). Note that the time scale is not 
logarithmic. 
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FIG. 4. Error in the electrostatic energy for the Na^ 2 cluster using the cubic cut-off method, as 
a function of its "control parameter" . Also shown is the time of calculation for each case (dashed 
line). Note that the time scale is not logarithmic. 
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